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1 . Introduction 


A 

Processes developed to date for calculating: the transonic poten- 
tial flow around airfoils and simple three-dimensional configurations 
depend for their evaluation on difference schemes which all require 
a rectangular coordinate system. 

For airfoils it is also important to set up a conformal pre- 
mapping to a circle where radii and concentric circles form a natural 
rectangular coordinate system. However, for the three-dimensional case 
such a process is hardly possible and besides, because of the after- 
differentiations would lead to quite complicated mathematical ex- 
pressions. 

The purpose of the present effort was to develop a method which 
operates exclusively directly in a physical plane so that nothing 
would prevent an extension to three dimensions. 

With certain limitations finite -element methods need not depend 
on orthogonal networks. However, they have a number of grave disad- 
vantages : 

- The solution depends on the element type. When the 
network is refined, methods with different elements 
converge to different results, and only in special 
cases do they lead to an exact solution. 

- The stream velocity is always underestimated. This is 
especially important for the transonic region where the 
magnitudes and the positions of compression shocks 
depend to a great degree on the velocity field. 

~ The proofs for convergence found in the literature are 
wrong . 
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- Althoup-h structural mechanics offer hundreds of element 
types, only one even approximately satisfies the accuracy 
requirements of progressive numerical aerodynamics. 

The last two items, in particular, create great difficulties in 
the development of the methods since one is forced to use computer- 
intensive cut-and-try methods. An extensive study has shown that, in 
a strictly mathematical sense the simple elements converge to an exact 
solution only when they are distributed across the flow field in an 
orthogonal combination. For then sm^h a method degenerates to a 
difference scheme. 

As always, the best aerodynamic element is the plane panel which, 
a priori, includes a non- trivial solution of the potential equation. 
However, it has not been possible to date to carry this over to the 
treatment of transonic flow problems. 

Thanks to a careful choice of the element and to the introduction 
of artificial viscosity, without which it would not at all be possible 
to achieve iterative convergence, the subject effort can at least be 
considered an engineering breakthrough of the finite-element method 
for transonic flow problems; however, the results shown in the litera- 
ture to date must be considered harmless, such as subsonic flow around 
circles, Joukowski airfoils and simple NACA airfoils v/hose field, in 
addition was often calculated in an approximately circular image plane, 
where the actual problems of the finite-element technique do not even 
come into play. 

2. Basic Equations ^ 

We shall base the mathematical formulation of our problem on the 
stationary, rotation-free, and isoenergetic flow of a completely ideal gas. 

There will be no sources if the continuity equation 

(gu)j< + (gw)jj .. 0 (1) 
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and there will be rotational freedom if 


MZ H w,< 


is satisfied I 


The energy carriesd alonf; by the fluid element is to be constant in 

the entire flow fields ? ? 2 

q<s • 




$ince we assumed irrotational flow, the changes of state are considered 
to be isentropici 

P " Po (4) 

For completeness sake we shall introduce the Mach number so that we 
can differentiate as to supersonic or subsonic flow: 


■ —rr 
^ i. t r* 


In (r; 


5 Po 

-X 7T^ T- 

Jo (1 » ~ 


3. Normalization 


Equations (3)* (^)i and (5) contain an unknown constant which 
can be easily removed if we define a dimensionless velocity q 
do cording -fco 


Then from equation (3) we get: 


4. ri2 

P08 2 ^ 




from (5) we f:eti 


( 7 ) 


. M? 

1 4 . m2 

substituting (4) into (6) gives j 

S“ 2o 

From (?) we can obtain the sonic velocity (M *= 1) to bes 



If we now understand that the uncancelled values represent the 
corresponding normalized values, then equations (1) and (2) remain 
unchanged. 

4, The Far-Field Solution / ^ 

4.1 Prandtl Transformation 

Independent of any further procedure the flow far downstream of 
the given airfoil can be estimated if equation (1) is made linear. 

Here we differentiate (1) tos 

The parallel onflow is disturbed only weakly so that in the following 
the vertical velocity component v/ can be neglected* 

♦ -0 ( 9 ) 

can be determined by after- differentiation from (8) 

q?o2 

’ 1 i.'”-' ?. 

1 - q 

2 

This can be written more simply if we replace q by equation (?) 

qjx » - r^2 j 
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( 10 ) 


At infinity the Kach numher will deviate only little from the incidsnt 
Kach number so that with (10) equation (9) simplifies to 

(1 - \^^?) Uj< ♦ W'2 « 0 , (11) 

Equation (2) permits the introduction of a potential according to 


so that (11) becomes 

'f’xx ♦ I.:: - 0 

From that we can form the Laplace equation is transformed with 


1 - H, 2 Z 

(14) 

This then produces! 

^xx ♦ 'rCf ■ 0 

(15) 


^l-*2 Integral potential representation 


If we substitute the expression 

X ♦ Z ♦ 


(16) 


into (15), then obviously the following equation must be solved! 

SV-x ♦ 7^5" 0 (17) 


To do this we use Green's theorem 

m 



VJe multiply (1?) by axi influence function e yet to be determined 
and integrate over (B) 

I 

JJ e<7xx+‘,Vi;)dx0r;^|J'oij-r,dr.- I ^ 7^)dxdj- 0 (10) 

(S) CKti^LC (S)“KLG 

■•A'" 


tjwg.* •A» 


»i i iiiiiirf 1-1^1^' 1 1 1 liii I iiii III III 



Sxohan'rln^^ e with 7 . we obtain analoi^ously 
(D)-l'^e CK*^LC (Q)-KLC 


The surface intef,ral on the ripht is replaced by (18) « /? 

§V«nds -JoifndR “ fl'f^Qxx <• oj^)dxd2 
CK-Le CK«»Ue (D)-KU& 

The contour pieces L and K make no contribution since the derivatives 
cancel each other over the forward and return path. 


Similarly the contribution of the outer border disappears at 
infinity since there ie pure parallel flow so that according to 
equation (16) the disturbance potential must disappear identically. 
Thus the followinf remains 

J^CnCis « Jo^nds ♦ f (ci^,,-»jGn)ds +J 

e e- c (n)-(i 


Tir'^ left integral can be evaluated immediately if a is small and if e 

is exclusively of function of the radius because then the integrand 

over the circle 6 is constant. 

, P.X 

B S '.J'p Oj. r d^ » ij’p Gj. r 2'!C 

>!l 

0 

if we now demand that 

Of r 2TT - 1 

then it follows that 


If we substitute this result into the integral equation and if we 
allow e to shrink toward 0, then we finally obtain 

fp “ In r - if (In r),,j dr. 

Where the integral must be taken over the edges of the airfoil. 
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^♦3 Blsplacemont Term 


We are seeking' an approKiraation solution for the first integral 
of (19) fL!r a slender airfoil. Then the following approximation 
hecomes valiiii 

^ - A 

iiri i)i: 

<l3 « dx 


or with (l^-) 


‘?D 


h'd - 4^ p <»« 


r (S In r d!< 


ZTcfl - j 


For disappearing supersonic velocities the linearized boundary 
condition 


/S 


is valid from which, by partial integration, we get 

<j?lj M Fb In r 1 - (5 z (In r)v dxl 

Zl p - 1. lo '» J 

For smooth closed profiles the first term drops out. Prom 
and equation (1^) we get 

(xp-x) dx 


‘fn 
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^ (xp-x) dx 

V . ■ — » 

(xp-x)^ + (zp~z)^ 


At infinity the profile irregularities at x, z become point-shaped 
so that in approximation we obtain 




(xp - x,n> 


2 F, (xp-Xrn)^ + (l~tv’) (zp-z,n)'’ J 


-^z clx (20) 


The singularity points x^, must he placed on the inside of the airfoil. 
The boundary integral becomes the airfoil surface. 




Mdv !gQrin 


lx 


Th® second integral in (19) can 139 evaluated as follows « 

** “ <Jo - 

fAf 


np 


< cp - s )^ + n ,(2 


- ■sVc -Ji?'' 


1 “=> 


- ^55^-i>d „to„ ^ 


- h § 


t?d^3' 


where d'Q' is the angle increment viewed from control point P, 



Partial integration* 


'Pr" ■ ^Vc ^ 


■k 


At a great distance from the airfoil the change of angle aisappears 
so that one can make the approximation 


fp- is § ■ il 'i'f - ^^5' 


V/ith (14) 


l|>p» r atari 




Xp - X 


m 


(21) 
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At the circulation slot the velocity is equal on both sides so that, 
surely, the following is trues 

o^> - o«(}» 

With the value for one obviously obtains 

‘{‘u •* " - r‘2-?: (’2) 

The potential jump thus is constant over the entire circulation slot. 


5 • Modification of the Sontinuitv Equation for T^gfeinsonic glow 


5*1 Dependence on Direction 


If at point P v/e calculate the velocity on which all gas-dynamic 
values depend according to our theory, then we approach the control 
point P from "both sides as we form the boundary value of the 
potential -differential quotient along the stream line s. 



The velocity thus always depends on a potential value ahead of and 
behind the control point, viewed in flov/ directiofi, regardless of 
how close one comes to the control point. 

For supersonic flow downstream of the control point signals 
are carried downstream at a faster rate than the expansion velocity 
of disturbances, flowing in the opposice direction, can carry them 
back to 'he control point. 

Thus the control point can only be influenced by events occurring 
upstream. Obviously for this case the described Cauchy potential 
derivation is not correct. On the contrary it is necessary that the 
velocity be generated at an auxiliary point H located upstream of the 
control point P. Here the distance H-P can be very small. 

Therefore we develop the following mathematical model.* /i2 

For subsonic flow we use the physical data needed for determining the 
potential as they are obtained for the control point itself. For 
supersonic flow we assign the physical data of a neighbor auxiliary 
point H lying upstream to the control point E-’itself. In this way the 



accuracy becomes correspondingly higher, as H comes closer to P. The 
contimiity eq.uation at P thus assumes the following forms depending 
on the type of flowi 

M < 1 " 0 
, H > t [(?v.)h]'..p t [<3 “>h]jp “ 0 


5.2 Artificial Viscosit?r 


For supersonic flow the described formulation can, using the 
values at the control point, be represented mathematically as a Taylor 
development which will first be demonstrated for (,*ii .* 


n («^u)p ♦ <Su)p3 Can - sp) 

formed as follows: 

‘S''>0 ■ '3;! q)p " (J|) jq ♦ J1 („q) 


The derivation is formed as follows: 


Since we move along a stream line in the small interval Sp-Sjj, its 
curvature plays no role so that disappears. 


The lowest arithmetical cost for a computer program is attained if, A 3 

instead of s we introduce the density g as independent variable: 

(gq)^ 4s - Cgq)g 4g » (q + gq^) 4^ « 

■ (q + ■§-) 40 » q (1 

aq J ‘ q§q •* 

With the aid of (10) we then obtain: 

(gq)j. Aa m q(l - 
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The expression in parentheses is negative for subsonic flow and dis- 
appears exactly for sonic velocity so that for differentiating the type 
of flow we can introduce t]ie switch function "max" 

(^u)|| " u max ^0,C 1 ~ 

Analogously* 

" j^p4- max [o,(l - - g)Jw 

Accordingly for transonic computer programs we need only replace the 

density by f i > 

2 ^ max (0,(1 - 

The additional methodology is not affected thereby. 

5*3 Viscosity Parameter /l4 

Since the length of the vector betv/een auxiliary point H and 
control point P is fixed by the discretionary dimension of the mesh 
network on which the computer program is based, the viscosity is 
generalized by the parameter G . 

J -r> 5 + e max [o,(l - 

Normally the value of e. will be close to 1. On the one hand, for 
a large number of meshes , the accuracy can be increased by the 
selection of a decreased viscosity parameter, but on the other hand 
there is danger that no iterative convergence is attained if the value 
of G. is too small . As shown by numerical experiments variations of C. 
can displace the compression shock within an interval of about 5f« 
of the airfoil height. 
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6t The Variation Principle 


6.1 Iv'eigrhted Residuum 


As 


Since we do not know, ahead of time, the solution of the continuity 
eq.uation (l) in conjunction with rotational freedom (12), (13), we have 
no choice except to introduce into it an approximating function with 
free parameters which must be chosen such that (1) is satisfied as well 
as possible. In this way the right side of equation (1) will generally 
not disappear exactly 

+ (gw);. • R 0 

The remainder R is called the residuum. 


There is a temptation to weight the residuum with a function in 
such a way that large residua are emphasized greatly while small 
residua are emphasized little and to demand that che integral thus 
formed disappear for the entire flow region. Such a weighting function 
is surely a deviation of the approximating potential distribution 
from the exact t 

j* f R ( (J) - dx dz o 0 


Here, however, is unknown. However one can, under the assumption 

that 9 deviates only little from the exact distribution, approximate 
the exact /p- values by the linear Taylor expansion* 


Introduce 


<j) => ^ + >, ( <J) - ) 

JJr dz a 0 


Little seems to be gained in this way since the unknown exact solution 
is still contained therein. However, the disappearance of the inte- 
gral can still be guaranteed if we demand that 

JJ R dx dz c 0 
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’A'lth (1), (12), (13) the followinf should then be true 

9 <J» ° 

Here, analogously to ( 18 ) partial integration can be used if we make 
the formal substitution* 

e — r» 

fx 

^ 9r|»S *)n - J*f g|?x f 2 dx d« • 0 

The first integral disappears since the mass-flux integral over the 
exterior boundary far downstream of the airfoil must be zero according 
to definition and since the fluid cannot penetrate the profile contour. 
If we interchange the differentiation in the second integral and write 
as in (12), (13) j then we reach the remarkably simple result 

[Jg(uu^ ♦ ww^ ) dx ds •• 0 

Thus the conto?.r boundary condition is implicitly contained therein. 

6.2 Minimum pressure integral 


Equation (23) can be further simplified. 

■ -rlh’ 


With ( 8 ) we get to 
1 C|^i|, dxdz » 



r* 





dxdz 


- 0 


The integrand is identified with ( 4 ) as a normalized pressures 


dx dz ■ 0 

This result can be confirmed in a simple manner from the principle 
of energy minima in closed systems. Imagine that a very large plate is 
dragged through a flow field where one side is subjected to the static 
pressure p of the flow while a constant equalization pressure is in 
effect on the back side. For the sake of simplicity let us assume that 


the plate is mounted vortically and that the drag direotion is hori 
zontal. !Then the work associated with the drag is 

A ■ J K fJx' •• bJJ(p - kontit') <Js dx 

The miniraum principle demands that A = A^^^, thus 

Jf(p - konst) dx dz « Min 


or 



0 


The result of this is tlu.t an exact minimum is never attained when an 
approximation is introduced. The pressure thus is overestimated while 
the velocity is underestimated. 


7 • Numerical Evaluation 


A 8 


?.l Choice of Element 

For the evaluation of (23) we require a stepwise approximation 
of the potential using free parameters which must be determined such 
that the value of the integral (23) disappears. This is customarily 
done by the method of finite <^lements . In deriving the minimum principle 
we started with the idea that the potential approximation does not 
deviate much from the exact solution. In contrast to the non-sensical 
convergence criteria often cited in the basic literature we state here 
precisely that: 

- The elements must be laid out in such a way that they 
always furnish the potential distribution. 

- For constant density the element must be source frees 

>• da - 0 

for all closed contours lying within the element including 
the boundary. This is the same as 

4^xx + " 0 


Therefore the density does not have to be considered 
because I as the size of the element becomes infinitely small, 
it can be placed as a constant in front of the throuyh-flow 
integral » 

- The element must have at least enoufih free parameters so 
that the validity of the continuity equation in inte/»;ral 
*<"orm can be assured, without contradiction, for any closed 
control contour. 

Except for the area panel W known in theoretical aero- 

dynamics there is no element which satisfies the above requirements. 
Since it has not been possible to date to incorporate the area panel 
into transonic computer programs, there is no other choice except to 
try to find, by trial methods, one of the known classical elements 
which comes closest to satisfy the engineering demands. 

a. The linear^triangle 

although with the equation 

^ ^ a ♦ bx ♦ cs 

it satisfies the continuity equation trivially, it does 
not satisfy the third convergence condition. One can 
easily see this, if somewhere in the flow field a velocity 
vector ^ and a single potential value <.)o is fixed. This 
is always possible by a proper choice of the constant of 
integration. 


In this way the potential values in 1 and ?, are known* 
We now add another triangle l-2-'3t 



If we set the mass flux at the separating line 1~2 equal /20 
on both sides, then we obtain the value i[> 3 . . If we com- 
plete the triangle configuration with the indicated side 3-0, 
then the mass flux condition via 2-3 ro longer furnishes 
the value 'J'o' • ^'he forced juncture-coupling of the tri- 
angles among each other leads to a redundancy. 

Figure 1 shows as an example a shock-free NI,R airfoil which 
was subsequently calculated with this element. As expected, 
the result is, as predicted, unusable, figure 2. 

The mutual coupling is cancelled only in the special case 
of right angle, isosceles triangles and the influence matrix 
for the potential values degenerates to that of the difference 
scheme. If the triangles are not isosceles, then 2 triangle 
sides must always coincide with a potential line and a stream 
line in order to achieve an exact convergence solution. 

Although Shen DJ and Habashi C'O report useful results with 
this element, they knowingly or unknowingly deceive the 
reader because the airfoil is first transformed into a 
near circular, conformal image plane where a nearly ortho- 
gonal network is constructed which additionally follows 
approximately the 4” resp. the lines. For non- ortho- 
gonal networks their methods are useless. 
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Brashears and Chan [5] work directly in the plane of flow 
and can p:ive only a rather challenfeahle interpretation 
of their weak results. Although Periaux and Perrier 0>3 
write a lot of theory, they produce no significant 
sensible results. 

Z he. au&dralic triaiig.^,e 

is based on the expression 

a ♦ bx * cr. ♦ dx?. ♦ ex^ ♦ 

Here the second convergence requirement is generally 
violated. The Laplace equation is satisfied only if ° ■ ~f 
This again is only true if the elements are arranged in 
a rectangular network. By an expansion to 6 parameters 
the juncture-coupling is only apparently cancelled com- 
pletely. In reality the juncture axes are now curved 
so that the redundancy is reduced, but still noticeable. 

For a method using this element the computer time is 
unjustifiably large. The unusable result shown in figure 3 
requires, e.g., 25 minutes on the central computer of 
an IBM 360/168. 


The f our:ippint_quadrangle. .el^me.nt 


is decidedly the best element. In a quadratic orthogonal 
combination the matrix degenerates to the difference scheme 
in the network rotated by 4-5° for that purpose. 


X I X X xr X X 

zliXixj X ! X i 

-'t'v ‘'i' '' I' 

y \ '>< I X J X I X 1 \ 
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In a right angle -quadrilattral comhination the Laplace 
equation is satisfied exactly everywhere. For an arbitrary 
element shape the Laplace equation is at least satisfied 
at th® Lhales circle via both radiating points. 

‘v’xx ♦ ^az " ® 


There is no redundancy. This element does not conform to 
the convergence criteria found in the text books, which 
incidentally are wrong} therefore it is seldom used, but 
satisfies the rctquirements listed here if the elements are 
distributed in such a way that the basic circles intersect 
the elements associated with them. 

d* The eighi-pointed quadraiigle_element 

is not as good as the four-pointed one since the matrix- 
diagonal elements for the corner points disappear exactly 
in the quadratic combination or otherwise are smaller than 
the neighboring elements so that the resulting system of 
equations for the potential values is conditioned extremely 
poorly. 

Hirsch [7] also reports relaxation factors smaller than 0.25 
so that we cannot consider this element. For a quadratic 
network division that author's method would fail for 
successive network refinement. 





ill 


Ab a basic ar«a 
use a sqtuare of 


for the approximation of a function f let us 
the following!* design i 

,iY<\ 


-1 




t 




-1 


f* 


let us say that the function f to be approximated is known at the 
corners 1 to Then an approximation is possible for the values 
of the function on the inside and at the ede:es of the square by 
usine the expression 


f • a ♦ >■ cr^ ♦ 


The coefficients are to be determined so that the polynomial 
takes on the exact values of the function at the corners! 

f^»a-b-c4d 
f2"a-b+c-d 
f3«a+b+c+d 
f/j^a + b- c- d 

Since the system of equations is nearly decoupled, the solution 
is simplified and after trivial 
followin^r result! , 

[fi 

♦ f2 

♦ f3 

+ u 


rearrangements leads to the 


(1 - 5 ) (1 -iq) 
(1 - g ) (1 ♦ ri) 
(1 4 g ) (1 

(1 4 g ) (1 - 1 ^)] 






'4/k 
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By alibre'-iatinr v,*e obtain the ^ollowinr /2^ 

f - rJ n Cl 
1 

whore is onllod the ‘f’orm function. 

7»3 S.te aMd£anrIc_el.ciiieiit 

The specialized results from the preceding; section can easily he 
carried ever to a fionoral qtuadranfle if this is transformed to 
the imaf^o aqtuare by means of the so-called isoparametoro. TI.xs 
is done by idontlfyinfr f successiYely as follows i 

•j) - I". 'Jl Ci (24) 

,1 

X - >- XI Cl (25) 

* • >. r.i Cl ( 26 ) 

1 



However, in the present minimization method the potential dis- 
tribution itself is not needed , but rather their derivations in 
the X- and s- directions s 


Tlx 


Now % and % cannot be inverted simply to functions of x and z 
so that we must first find the distortion determinant as a function 
of ij and . Here v/o v/rite the total differential ass 

ill 

dj » dx ♦ dz 
di^ " Tlx dx ♦ Tl-;; dz 
dx - x^ dg + x-q^ d'fi 
da . .zg dg + d-q^ 
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± 1 ^ 






If WG insert the upper two equations into the expressions for dx 
and da, then we obtain 

dx - <xg 5;{ ♦ xn^'n,<)c5H ♦ ♦ X5 1ia>d» 

. dz • (z|i ♦ <=qHz ♦ 5’f, ^»)dz 


dx must be independent of da so that by exchanging coefficients 
we obtain four equations! 

Ijx ♦ lx • 1 

Iz ♦ Xg • 0 

®>l ♦ zg Ta - 1 

“Sx ♦ Z7^ lx • 0 


Solving for the unknowns ?x» ?a» lx* la 

Sx ■ 


5?/- 

u 


il?L 
" u 


(27) 

( 20 ) 


'lx ■ 
Iz " 


’’*1 

D 



(29) 

(30) 


Wjl^A 0 ■ xg " x^ Z'g (31) 

Thus all operations for the formation of the functional relations 
have been prepared. 


7.4 jEnergy_of Ihe element /2g 

The term element energy is defined as the value of the integral (23) 
over the element. As a pivotal point we choose the corner 1. In 
order to keep the computer cost as low as possible we assume that 
the density is constant as distributed over the element. Using 
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the chain rule of differentiation we obtain, with equations (27) 
to (3i) the following expression for the element energy* 


I 



- Gig 



Tt f, (C 






“g’J J “ i5~ 


(32) 


7.5 Numerical integration 

It is extremely difficult to determine the value of the integral 
(23) in an elementary manner. Therefore one tries to find at 
least a good approximation value using Gaussian integration. For 
the one-dimensional integral we first form the expression 

1 

f f dg - a f I'ii) ♦ b f ( 1 ^ 2 ^ 


This formula should be exact for the four test functions 

fn - g" Cn - 0, i, 2, 3) 


OF 


Then the following is true 


n ■ 0 


n ■ 1 


n ■ 2 


n » 3 


1 

J dj; •• 2 «• ,1 b 


1 

J u 0 « ajji + bj?2 

-1 





The thus resulting system of equations for the four unknowns 
'Ell Ji 2 can he solved in a simple way and leads to the 
following result k i ‘ 


El " 
^2 " 


1‘ 



1 

iT 


Thus for the general integral the following is truei 




The double integral can be evaluated as an approximation by means 
of 2-fold Gaussian quadratics! 


-1 -1 
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1 1 


1 1 


1 1 


f(- -i,— '-) + f(-4— -i-) + f(-i- _i_)^ 

^ f3" i3 iTii' 


7*6 Stiffness matrix 
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7.6.1 Pield-£.olnt collocation 


Until now the special corner 1 of an element was chosen as the 
starting point for the element energy. It is understood that in a 
global element combination the four elements coming together at a 
point furnish a contribution to the integral (23). Here (32) is con- 
sisddred to be unchanged if the element corner numbering system is 
retained locally by a corresponding rotation of the coordinates* 



Since we use column relaxation in the computer program, instead of the 
entire stiffness matrix we need only build up the tri -diagonal matrix 
for a selected line. The contributions of discretionary points which 
are not located on this line, are placed on the right-hand side of the 
system of equations* 



o matrix contribution 
X right-hand side 
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7.6.2 Boundary_j)oint_coliocation 


At the extreme boundary of the computer network the potential 
values are assumed for the far-field solution ( 16 ), ( 20 ), ( 21 ) 



o matrix contribution 
^ right-hand side 


At the boundary of the airfoil obviously only two instead of four 
elements contribute to the formation of the column matrix. 



7 . 6.3 Girculat ion-slpt_c ollocati, 


on 


To determine the circulation r for equation (21) the potential 
values the potential values are calculated directly aoccrding to ( 22 ) 
after each field relaxation, at the upper and lower surface of the 
trailing edge of the airfoil. 


original page 13 

OF POOR QUALITY 
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For the slot itself we make use of ( 22 ) and, e.g., replace <f»o with 
(Vi^j + . Only then can the derivation be carried out according to 

(23) from the potenti.al values of the lower bank. As in field 
collocation, the contributions of four elements each constitute a 
pivotal value. The contributions of the element energy which were 
multiplied by p are transposed to the right-hand side of the tri- 
diagonal matrix! 



matrix 

° contribution 

X right-hand 
side 


8. Results 


Z2I 


WXR lifting Airfoil 


As already mentioned, not all elements give the same 
results. The disadvantages of the linear triangle have 
already been discussed. Thus it is only necessary to 
look at figures 1 and 2 to discover that this element 
cannot be used. 
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For similar reasons the quadratic triangle also fares 
poorly in the numerical comparison, figure 3* 

As an experiment the quadrangle element was tested to see 
if it could satisfy the continuity equation in a conser- 
vative manner. As shown in figure 4 the pressure is greatly 
underestimated. A similar procedure is suggested by 
Jameson and Caughey DOtV/ho, however, seem to obtain 
usable results. Through refinement of the network the 
mentioned effect is enhanced, figure 5* Only the 
quadrangle element in conjunction with the principle of 
variations produces a satisfactory result, figure 6. 

- NAG A 0012 

To accelerate convergence the following steps were taken s 

- direction of relaxation •= main stream direction 

- excess relaxation of IQOfo for 

~ intensive subdivision of the network 

/32 

A series of networks and results is shown in figures 7 "fco l6. 

For a transition to whatever network follows the potential 
distribution is interpolated linearly. The mesh buidup 
is controlled by the parameters and, of course, is done 
automatically. 

The end result, figure l6, is not satisfactory since the 
viscosity parameter was chosen too small by QOfo. For 100^ 
the strength of the shock becomes realistic, figure I 7 . 


- NAG A 64A410 


®his result is in p-ood aitreement with those in 
standard publications, figure 18. 

- Whitcomb airfoil , figure 19 

The profile generated automatically according to [ 9 ] 
was to show to what degree strong "rear loading" is 
reproduced. Here we must point out the necessity for 
Gaussian integration to obtain the element energy. 

The version of the program using trapezoidal inte- 
gration fails for this example. 

- Korn airfoil 

Figure 20 points out the importance of the viscosity 
parameter. Although the solution does converge, the 
viscosity parameter for this highly sensitive profile 
must be set at 200 ?^ in order to dampen out the waviness 
in the supersonic region, figure 21 ; this is an indi- 733 

cation that the method performs too "stably" so that, 
by increased viscosity, the correct root can be found 
from the many possibilities of the differential 
equations . 


As the result of the steps taken to accelerate convergence 
the computer times are comparatively small. After 20 iterations 
on the last network the pressure distribution was entirely frozen. 

This produced computer times for the central unit of an IBM 360/168 

1 minute for the 62 . 16 network 
4 minutes for the 124 . 32 network 
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9. Summary 


sTlith the introduction of a special method of finite elements in 
conjunction with the mathematical model of artificial viscosity it 
has been possible to produce a usable computer method for transonic 
flow around airfoils which does not require an orthogonal computer 
network so that the method can, in a simple manner, Je adapted to 
three-dimensional disturbance problems. Such a program is in a 
test stage at the present time. 
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